Supplement to ‘A comparative study of network-based pathway enrichment analysis methods’

Test Design

As described in Section 3 of the main paper, three different pathway dysregulation mechanisms were used, whose characteristics are summarized in Table 1. Specifically, for all topology designs, we set the detection call (DC) to be 10%, meaning that 10% of the genes/metabolites inside the pathway are differentially expressed with a mean difference varying from 0.1 to 0.5. Note that the magnitude of the mean signal is expressed relative to the unit variance of each gene/metabolite. Further, at each mean level, we examine the results of enrichment analysis methods with and without permuting the sample labels. When the sample labels are permuted, data from the two experimental conditions share a similar correlation structure. As a result, the significance of each pathway is primarily driven by the differential expression of its members. On the other hand, when the sample labels are fixed to be the same as the ones in the study, it is expected that the two experimental conditions may have different correlation structures. Therefore, enrichment of each pathway could be due to differential expression and/or differences induced due to changes in their interactions.

Table 1: Design of simulation experiments at 10% DC

Topology.design	Type	Permutation	Mean
Community	1	no	0.1
			0.2
			0.3
			0.4
			0.5
	2	yes	0.1
			0.2
			0.3
			0.4
			0.5
Neighborhood	3	no	0.1
			0.2
			0.3
			0.4
			0.5
	4	yes	0.1
			0.2
			0.3
			0.4
			0.5
Betweenness	5	no	0.1
			0.2
			0.3
			0.4
			0.5
	6	yes	0.1
			0.2
			0.3
			0.4
			0.5

Results

In this section, we present numerical results with different pathway dysregulation designs for the two cancer studies and for the metabolomic animal study, in addition to the ones shown in Section 4 of the main paper. Due to space limitation, we show only the average powers across multiple pathways as done in the main paper. Detailed results regarding the empirical powers for each pathway are available on GitHub.

Analysis of breast cancer data

We start with comparisons in the TCGA breast cancer study with community and neighborhood dysregulation designs, together with more pathway level networks.

Empirical powers

The average powers for pathways grouped by DC and pathway size with varying mean signals are shown in Figure 1 using the original sample labels from the TCGA breast cancer study and Figure 2 using shuffled sample labels under the community dysregulation design.

Figure 1: Average powers of multiple pathways, grouped by DC and pathway sizes, using sample labels from the original study and based on community dysregulation, in the TCGA breast cancer study (Koboldt et al. 2012). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. No method is consistently better than its competitors across different mean levels when DC is small (top two panels), whereas PathNet and CAMERA show minor advantages when DC is large (bottom two panels).

Figure 2: Average powers of multiple pathways, grouped by DC and pathway sizes, using shuffled sample labels and based on community dysregulation, in the TCGA breast cancer study (Koboldt et al. 2012). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. Overall, DEGraph outperforms the others, with CAMERA exhibiting the worst performance.

The average powers for pathways grouped by DC and pathway size with varying mean signals are shown in Figure 3 using the original sample labels from the TCGA breast cancer study and Figure 4 using shuffled sample labels under the neighborhood dysregulation design.

Figure 3: Average powers of multiple pathways, grouped by DC and pathway sizes, using sample labels from the original study and based on neighborhood dysregulation, in the TCGA breast cancer study (Koboldt et al. 2012). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. No method is consistently better than its competitors across different mean levels when DC is small (top two panels), whereas PathNet and CAMERA show minor advantages when DC is large (bottom two panels).

Figure 4: Average powers of multiple pathways, grouped by DC and pathway sizes, using shuffled sample labels and based on neighborhood dysregulation, in the TCGA breast cancer study (Koboldt et al. 2012). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. Overall, DEGraph outperforms the others, with CAMERA exhibiting the worst performance.

Overall, we observe that the performance of PathNet is better when using fixed sample labels, whereas DEGraph performs better when the sample labels are permuted because the underlying assumption of same network structure across the two conditions under consideration is satisfied. In contrast, CAMERA exhibits poor performance when the sample labels are permuted, since the gene-gene correlations between the two experimental conditions are approximately the same. This likely implies that the way CAMERA characterizes the gene level differential expression is suboptimal compared to the competitive test in PathNet. In all settings, the performance of NetGSA is reasonably reliable and is comparable to DEGraph and PathNet in most settings.

Network of enriched pathways

In Section 4 of the main paper, the network of enriched pathways (defined as those with empirical powers at least 0.8) under the betweenness dysregulation at mean signal 0.5 was presented to highlight that NetGSA, DEGraph and PathNet have high concordance for enrichment analysis of genetic pathways. For completeness, we provide the networks of enriched pathways under the community and neighborhood dysregulation designs at mean signal 0.5. Similar networks at other mean signal levels can also be constructed.

Figure 5 and Figure 6 show, respectively, the network of enriched pathways without and with sample permutations under the community dysregulation design. The two figures resemble each other and confirm that the three methods identify a similar subset of enriched pathways regardless of sample permutations.

Figure 5: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on community dysregulation and fixed sample labels for the TCGA breast cancer study (Koboldt et al. 2012). Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared, out of the total unique genes in the two pathways.

Figure 6: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on community dysregulation and shuffled sample labels for the TCGA breast cancer study (Koboldt et al. 2012). Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared, out of the total unique genes in the two pathways.

Similarly, Figure 9 and Figure 10 show, respectively, the network of enriched pathways without and with sample permutations under the neighborhood dysregulation design. The two figures resemble each other and confirm that the three methods identify a similar subset of enriched pathways regardless of sample permutations.

Figure 7: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on neighborhood dysregulation and fixed sample labels for the TCGA breast cancer study (Koboldt et al. 2012). Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared, out of the total unique genes in the two pathways.

Figure 8: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on neighborhood dysregulation and shuffled sample labels for the TCGA breast cancer study (Koboldt et al. 2012). Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared, out of the total unique genes in the two pathways.

In summary, under the community and neighborhood design, we also observe high concordance among NetGSA, DEGraph and PathNet when studying genetic pathway enrichment from gene expression data.

Analysis of prostate cancer data

Type I errors

Table 2 shows the average type I errors of multiple pathways, grouped by pathway sizes, in the TCGA prostate cancer study (Abeshouse et al., 2015). Again, we notice the similar pattern as in the TCGA breast cancer study: topologyGSA has highly inflated type I errors, whereas CAMERA and DEGraph both have very tight control of type I errors. NetGSA and PathNet have comparable performances.

Table 2: Average type I errors over multiple pathways, grouped by pathway sizes, in the TCGA prostate cancer study (Abeshouse et al. 2015)

	Pathway size
Method	<=75	>75
NetGSA	0.024	0.061
DEGraph	0.001	0.001
topologyGSA	0.283	0.570
SPIA	0.000	0.000
Pathway-Express	0.000	0.000
CAMERA	0.007	0.006
PathNet	0.034	0.054

Empirical powers

Figure 9 and Figure 10 depict the average powers for multiple pathways in the TCGA prostate cancer study under the betweenness dysregulation design.

Figure 9: Average powers of multiple pathways, grouped by DC and pathway sizes, using sample labels from the original study and based on betweenness dysregulation, in the TCGA prostate cancer study (Abeshouse et al. 2015). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. No method is consistently better than its competitors across different mean levels when DC is small (top two panels), whereas PathNet and CAMERA show minor advantages when DC is large (bottom two panels).

Figure 10: Average powers of multiple pathways, grouped by DC and pathway sizes, using shuffled sample labels and based on betweenness dysregulation, in the TCGA prostate cancer study (Abeshouse et al. 2015). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. Overall, DEGraph outperforms the others, with CAMERA exhibiting the worst performance.

Figure 11 and Figure 12 depict the average powers for multiple pathways in the TCGA prostate cancer study under the community dysregulation design.

Figure 11: Average powers of multiple pathways, grouped by DC and pathway sizes, using sample labels from the original study and based on community dysregulation, in the TCGA prostate cancer study (Abeshouse et al. 2015). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. No method is consistently better than its competitors across different mean levels when DC is small (top two panels), whereas PathNet and CAMERA show minor advantages when DC is large (bottom two panels).

Figure 12: Average powers of multiple pathways, grouped by DC and pathway sizes, using shuffled sample labels and based on community dysregulation, in the TCGA prostate cancer study (Abeshouse et al. 2015). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. Overall, DEGraph outperforms the others, with CAMERA exhibiting the worst performance.

Finally, Figure 13 and Figure 14 depict the average powers for multiple pathways in the TCGA prostate cancer study under the neighborhood dysregulation design.

Figure 13: Average powers of multiple pathways, grouped by DC and pathway sizes, using sample labels from the original study and based on neighborhood dysregulation, in the TCGA prostate cancer study (Abeshouse et al. 2015). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. No method is consistently better than its competitors across different mean levels when DC is small (top two panels), whereas PathNet and CAMERA show minor advantages when DC is large (bottom two panels).

Figure 14: Average powers of multiple pathways, grouped by DC and pathway sizes, using shuffled sample labels and based on neighborhood dysregulation, in the TCGA prostate cancer study (Abeshouse et al. 2015). The x-axis shows the magnitude of mean difference added to the a ected genes. In general, power increases as the mean difference increases. Overall, DEGraph outperforms the others, with CAMERA exhibiting the worst performance.

Overall, the patterns are similar to those observed in the breast cancer study. NetGSA shows robust performance throughout. In comparison, PathNet performs better when the sample labels are fixed, whereas DEGraph performs better when sample labels are permuted. SPIA and Pathway-Express (PE) have similar performances and only work for larger mean signals.

Network of enriched pathways

Lastly, the networks of enriched pathways in the TCGA prostate cancer study under the betweenness design are presented in Figure 15 with fixed sample labels and Figure 16 with shuffled sample labels. The two figures resemble each other and confirm that the three methods identify a similar subset of enriched pathways regardless of sample permutations.

Figure 15: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on betweenness dysregulation for the TCGA prostate cancer study (Abeshouse et al. 2015), using fixed sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared out of the total unique genes in the two pathways.

Figure 16: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on betweenness dysregulation for the TCGA prostate cancer study (Abeshouse et al. 2015), using shuffled sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared out of the total unique genes in the two pathways.

The networks of enriched pathways in the TCGA prostate cancer study under the community design are presented in Figure 17 with fixed sample labels and Figure 18 with shuffled sample labels. The two figures resemble each other and confirm that the three methods identify a similar subset of enriched pathways regardless of sample permutations.

Figure 17: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on community dysregulation for the TCGA prostate cancer study (Abeshouse et al. 2015), using fixed sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared out of the total unique genes in the two pathways.

Figure 18: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on community dysregulation for the TCGA prostate cancer study (Abeshouse et al. 2015), using shuffled sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared out of the total unique genes in the two pathways.

The networks of enriched pathways in the TCGA prostate cancer study under the neighborhood design are presented in Figure 19 with fixed sample labels and Figure 20 with shuffled sample labels. The two figures resemble each other and confirm that the three methods identify a similar subset of enriched pathways regardless of sample permutations.

Figure 19: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on neighborhood dysregulation for the TCGA prostate cancer study (Abeshouse et al. 2015), using fixed sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared out of the total unique genes in the two pathways.

Figure 20: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 based on neighborhood dysregulation for the TCGA prostate cancer study (Abeshouse et al. 2015), using shuffled sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether PathNet, DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of genes shared out of the total unique genes in the two pathways.

In summary, we see the same concordance among NetGSA, DEGraph and PathNet in detecting enriched pathways, as observed in the TCGA breast cancer study.

Analysis of metabolomic data

The network of enriched metabolic pathways with different mean signals and sample labels are shown, respectively, in Figure 21 and Figure 22. Consistent with the power curves shown in the main paper, NetGSA reports more enriched pathways when using fixed sample labels, whereas DEGraph yields more enriched pathways when using permuted sample labels.

Figure 21: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 for the metabolomic study (Fahrmann et al. 2015), using fixed sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of metabolites shared out of the total unique metabolites in the two pathways.

Figure 22: The network of enriched pathways (empirical power at least 0.8) with mean signal 0.5 for the metabolomic study (Fahrmann et al. 2015), using permuted sample labels. Each node represents a pathway, with node size proportional to the size of the pathway. Nodes are colored based on whether DEGraph and NetGSA identify the pathway as enriched. An edge is drawn between two pathways if they have one or more genes in common. Only the top 5% of all edges are visualized. Width of the edge indicates the proportion of metabolites shared out of the total unique metabolites in the two pathways.

Analysis of synthetic data

Lastly, we present a simulation example based on synthetic data to illustrate the advantage of NetGSA in detecting both changes in mean expression levels and/or changes in network structure. Figure 23 shows the network and subnetworks being compared.

Figure 23: The network and subnetwork topology under the null (left) and alternative (right). Dashed lines represent edges that are present in only one condition. Nodes in square are associated with mean changes.

As indicated by node shapes, subnetworks 1 and 6 have equal means under the two conditions, whereas subnetworks 3 and 8, 4 and 5, 2 and 7 have, respectively, 20%, 40% and 60% nodes with differential means. Further, some subnetworks also have changes in network topology. In particular, the topologies of subnetworks 4, 6, 7 and 8 remain the same under both conditions, whereas those of subnetworks 1, 2, 3, and 5 are completely rewired. Such dramatic changes in network topology are often expected in practice. In this case, DEGraph is not recommended because it assumes the network structures under the two conditions are the same, when in fact they are not. Whichever network provided for DEGraph will misspecify the graph topology for at least one condition. In contrast, NetGSA can leverage both changes in mean expressions and changes in network structures.

The empirical powers for each subnetwork averaged in 100 replications are shown in Table 3. Here the results for NetGSA were obtained assuming the network under each condition is known a priori. The network provided for DEGraph is the network under the null. As expected, NetGSA identifies 5 (0.99) and 3 (0.96) as the most significantly enriched pathways due to both changes in mean and in network topology, followed by pathway 7 (0.94) and 2 (0.89). An interesting observation is that NetGSA assigns slightly higher power for pathway 7 compared to pathway 2. This is because some of the DE nodes in pathway 7 are in hub positions. In comparison, DEGraph correctly identified pathway 2, 3, and 5 as the enriched ones, but missed pathway 7.

It is important to note that it is unfair to compare NetGSA and DEGraph using this synthetic data example when the network for DEGraph has to be misspecified. Nonetheless, this example illustrates how NetGSA can be used to simultaneously detect changes in mean and network structures.

Table 3: Empirical powers averaged in 100 replications with 50 observations per condition

Method	set.1	set.2	set.3	set.4	set.5	set.6	set.7	set.8
NetGSA	0.08	0.89	0.96	0.14	0.99	0.02	0.94	0.03
DEGraph	0.18	1.00	1.00	0.49	1.00	0.06	0.62	0.31
true power	0.12	0.93	0.98	0.11	0.99	0.05	0.95	0.10

R Code

For reproducibility, R code used in this comparative study is available on GitHub.

References

Abeshouse, Adam, Jaeil Ahn, Rehan Akbani, Adrian Ally, Samirkumar Amin, Christopher D Andry, Matti Annala, et al. 2015. “The Molecular Taxonomy of Primary Prostate Cancer.” Cell 163 (4). Elsevier: 1011–25.

Fahrmann, Johannes, Dmitry Grapov, Jun Yang, Bruce Hammock, Oliver Fiehn, Graeme I Bell, and Manami Hara. 2015. “Systemic Alterations in the Metabolome of Diabetic Nod Mice Delineate Increased Oxidative Stress Accompanied by Reduced Inflammation and Hypertriglyceridemia.” American Journal of Physiology - Endocrinology and Metabolism 308 (11). American Physiological Society: E978–E989.

Koboldt, Daniel C, Robert S Fulton, Michael D McLellan, Heather Schmidt, Joelle Kalicki-Veizer, Joshua F McMichael, Lucinda L Fulton, et al. 2012. “Comprehensive Molecular Portraits of Human Breast Tumours.” Nature 490 (7418). Nature Publishing Group: 61–70.